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I.  INTRODUCTION 


In  Part  1  of  this  article  we  studied  the  steady-state, 
one-dimensional  solutions  of  an  ionizing  shock  in  Argon 
using  a  detailed  collisional-radiative  model.  In  this  sec¬ 
ond  part,  we  focus  on  the  unsteady  and  two-dimensional 
features  of  the  problem.  While  there  have  been  nu¬ 
merous  studies  of  ionizing  shock  propagation  in  noble 
gases,  it  was  not  until  the  shock  tube  experiments  con¬ 
ducted  at  the  Institute  for  Aerospace  Studies,  Univer¬ 
sity  of  Toronto  (UTIAS),1  that  it  was  shown  that  under 
certain  conditions,  a  definite  dynamical  behavior  could 
be  associated  with  a  strong  ionizing  shock  in  argon.  It 
was  discovered  that  the  translational  shock  front  prop¬ 
agating  in  pure  argon  “develops  sinusoidal  instabilities 
which  affect  the  entire  shock  structure  including  the  ion¬ 
ization  relaxation  region,  the  electron  cascade  front,  and 
the  final  quasi-equilibrium  state.”  Deviations  from  ideal, 
one-dimensional  shock  solutions  are  fequently  found  in 
practice,  due  to  finite-size  effects  in  experimental  de¬ 
vices.  Thus,  shock  curvature  and  boundary-layer  ef¬ 
fects  in  ionizing  shocks  can  be  expected.  For  example, 
Bristow2  postulated  the  presence  of  ”weak  forward  mov¬ 
ing  oblique  shocks”  while  Brimelow3  showed  that  in  pure 
argon,  premature  ionization  occured  close  to  the  wall 
(l-2mm),  which  can  be  easily  explained  from  the  fric¬ 
tional  heating  in  the  developing  boundary  layer.  How¬ 
ever,  Bristow  also  noticed  a  multiply  curved  shock  front 
with  what  appeared  to  be  transverse  waves,  an  effect  he 
conjectured  could  be  caused  by  amplification  of  acoustic 
perturbations  through  thermal  non-equilibrium  effects. 
Buchanan  performed  time-resolved  interferometry  as  well 
as  Schlieren  imagery  and  noticed  transverse-moving  and 
colliding  shocks.  Thus,  he  obtained  images  with  cross- 
hatching  features  which  resembled  the  pattern  left  by 
a  multi-headed  gas-detonation  wave  on  a  soot-covered 
wall.  He  also  computed  the  cell  size  from  the  number  of 


half-sine  waves  that  were  present  across  the  translational 
front.  Houwing  et  al.A  carried  out  shock-tube  experi¬ 
ments  under  similar  conditions  to  the  UTIAS  cases  and 
observed  similar  instabilities,  albeit  not  as  a  dramatic. 
However  their  shock  tube  had  a  much  smaller  cross  sec¬ 
tion  with  a  diameter  of  5cm,  compared  to  the  10x18  cm 
shock  tube  at  UTIAS.  The  presence  of  impurities  can  also 
affect  the  nature  of  the  relaxation  region.  Bristow  showed 
that  adding  a  small  amount  (0.4%  partial  pressure)  of  H2 
was  sufficient  to  make  the  shock  planar.  Therefore,  the 
absence  of  instabilities  in  some  experimental  results  may 
very  well  be  a  result  of  length  scales  and  experimental 
conditions,  notably  gas  purity. 

Drawing  from  the  summarizing  work  of  Glass  &  Liu1, 
Cambier5  studied  the  transient  effects  of  such  shocks  in 
a  series  of  numerical  experiments,  verifying  the  concept 
of  a  “resonant  cavity”  phenomena  proposed  by  Bristow. 
Modeling  the  same  conditions  as  the  UTIAS  experiments 
he  was  able  to  obtain  self-sustaining  oscillations  of  the 
shock  front  that  were  similar  to  the  galloping  instabili¬ 
ties  of  detonation  waves.  It  was  shown  that  the  insta- 
bilites  were  sensitive  to  the  collisional-radiative  kinetics 
and  that  the  instabilites  would  be  damped  if  the  exci¬ 
tation  and  ionization  cross  sections  fell  outside  a  given 
range.  He  also  verified  the  effect  of  molecular  impu¬ 
rities  in  the  dampening  of  instabilities.  However,  that 
work  was  performed  only  in  one  dimension  and  with  a 
very  simplified  collisional-radiative  model;  the  relation¬ 
ship  with  the  cell  pattern  experimentally  observed  could 
only  be  conjectured,  although  Cambier  was  able  to  pre¬ 
dict  the  ionization  cell  size  to  a  satisfactory  degree.  We 
extend  here  this  study  to  two  dimensions  in  order  to 
clearly  demonstrate  the  mechanism  responsible  for  in¬ 
stabilities  in  ionizing  shocks  in  pure  noble  gases.  Since 
instabilities  were  observed  in  both  two  and  three  dimen¬ 
sions,  it  is  sufficient  to  perform  these  computations  in 
(planar)  2D — extension  to  3D  unsteady  calculations  is 
contemplated  for  the  near  future. 
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II.  GOVERNING  EQUATIONS 
A.  Plasma  thermodynamics 

We  describe  the  plasma  as  a  neutral  single- fluid,  two- 
temperature  flow.  The  individual  electronic  states  of  neu¬ 
tral  argon  are  separately  convected  as  different  species, 
allowing  a  non-Boltzmann  description  of  the  Atomic 
State  Distribution  Function  (ASDF).  Charge  neutrality 
is  guaranteed  by  the  small  Debye  length  scale  and  the 
absence  of  applied  electric  fields,  while  a  single  fluid  ve¬ 
locity  is  a  consequence  of  the  large  effective  frequency 
for  collisional  momentum  exchange  between  the  various 
plasma  components,  for  the  typical  densities  of  interest 
here.  Presently,  we  neglect  viscous  and  species  diffusion 
effects.  In  a  finite- volume  framework,  the  hyperbolic  sys¬ 
tem  of  equations  has  the  following  form: 

~r  [  Q  dV  +  (f  h  ■  F  dA  =  f  Cl  dV  (1) 
dt  Jv  J A  Jv(t ) 


i.e.  the  partial  derivatives  are  evaluated  for  all  other  con¬ 
served  variables  being  constant  (see  Appendix  A).  The 
derivatives  of  the  pressure  are  thus  obtained  and  using 
the  notation  PQi  =  J^,  we  obtain  the  speed  of  sound: 

a2  =  £  JPp.  +  (£  +  ~p  -  y )  pe  +  SePSe  (4) 

The  hyperbolic  system  of  the  Euler  equations  has  the 
set  of  eigenvalues  vn,  vn  db  a  and  Q  and  F  can  be  decom¬ 
posed  according  to  the  left  and  right  eigenvectors  (see 
A).  Thus  the  Jacobian  is: 

a=(%)=ral  (5» 

The  term  Cl  contains  the  source  terms,  both  collisional 
and  radiative,  and  is  described  in  complete  detail  in  ar¬ 
ticle  1. 

B.  Numerical  Scheme 


where  V  is  the  volume  of  an  elementary  fluid  ’’cell”,  A 
is  the  bounding  area  and  h  its  normal.  Q  is  the  vec¬ 
tor  of  conserved  variables,  F  is  the  flux  and  Cl  contains 
the  source  terms  due  to  the  kinetics  of  the  collisional- 
radiative  processes.  The  conserved  variables  and  the  nor¬ 
mal  flux  can  be  written  as: 


The  system  is  solved  using  the  operator-splitting  ap¬ 
proach:  at  each  time  step  the  convective  terms  are  solved 
independently  of  the  source  terms  and  the  conserved  vari¬ 
ables  are  updated  for  each  cell  via 

Qn+1  =  Qn  +  AQcqnv  +  AQcr  (6) 
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with 

E  =  Eh  +  Ee  +  -pF2  (3) 

the  total  energy  density  of  the  plasma  and  H  =  E  +  p 
the  total  enthalpy  density  and  p  =  Ph  +pe  is  the  total  gas 
pressure.  We  define  Se  =  pse,  with  se  =  Pe/p7e  as  an 
entropy-like  variable,  used  instead  of  the  thermal  energy 
of  the  electron  in  order  to  maintain  a  strictly  hyperbolic 
formulation6.  This  is  necessary  for  accurate  determina¬ 
tion  of  the  electron  temperature  across  the  shock  discon¬ 
tinuity.  Note  also  that  we  use  the  total  mass  density  in 
the  definition  of  Se  instead  of  pe,  to  avoid  singular  values 
in  the  limit  of  vanishing  ionization  fraction.  Total  dif¬ 
ferentials  of  the  conserved  variables  can  be  obtained  (see 
Part  1): 


Because  of  the  short-time  scales  involved,  the  collisional- 
radiative  source  terms  are  updated  implicitly, 

AQcr  =  A t(I  -  JAt)^Cl(Qn)  (7) 

requiring  computation  of  the  Jacobian,  J  =  dCl/dQ ,  and 
the  use  of  a  Gauss- Jordan  elimination  step.  The  con¬ 
vective  terms  are  updated  by  the  second-order  Adams- 
Bashforth  integration  scheme,  which  uses  the  fluxes  F(Q) 
computed  at  a  time-interpolated  value  of  the  conserved 
variables,  using  current  and  previous  solutions: 

AQconv  =  -y  £  *A.(§Q"  ~  \Qn~X)  ^  (8) 


the  summation  being  over  all  faces  of  a  computational 
cell,  with  normal  vectors  ns  directed  outwards.  An  ad¬ 
vantage  of  the  Adams-Bashforth  scheme  is  the  fact  that 
the  fluxes  need  to  be  evaluated  once,  as  opposed  to  a 
multi-step  method  like  Runge-Kutta;  this  also  reduces 
the  overhead  of  message-passing  in  distributed-memory 
architectures. 

The  term  F  in  Eq.  8  denotes  the  numerical  flux  which 
we  compute  using  the  hlle  Riemann  solver7, 
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The  upwind  wave  speeds  b  are  given  by 


MQ)  =  E 


b+  =  max(0,  v„,R  +  aR,  vn  +  a)  (10) 

b~  =  min(0,  v„)L  -  aL,  vn  -  a)  (11) 
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where  the  tilde  symbol  (7)  is  used  to  denote  the  Roe- 
average  value  at  the  cell  interface.  High-order  spatial 
resolution  is  achieved  via  parabolic  interpolation  of  the 
left  and  right  states  at  the  cell  interfaces  via 

Ql  =  -(2Qj_i  +  5  Qj  —  Qj+ 1)  (12) 

making  the  scheme  3rd-order  accurate8.  We  use  the  bar 
(7)  symbol  to  differentiate  cell-average  quantities  from 
point  values.  Strong  nonlinear  waves  require  limiting, 
which  modifies  the  left  and  right  states  according  to 

Ql  <-  median(QL,  Qj ,  QMP)  (13) 

with 

QMP  =  Qj  +  minmod[QJ+i  -  Qj ,  a(Qj  -  Qj- 1)]  (14) 

being  the  monotonicity-preserving  (MP)  limit9.  Here,  we 
choose  the  parameter  a  to  be  2.  The  right  state  is  easily 
found  from  symmetry.  Due  to  the  non-linearity  of  the 
limiter,  the  reconstruction  is  performed  on  the  charac¬ 
teristic  variables,  requiring  a  full  diagonalization  of  the 
governing  equations. 


III.  ID  RESULTS 
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FIG.  1.  Numerical  setup  for  unsteady  ID  shock  tube  simula¬ 
tions. 


FIG.  2.  Instantaneous  Mach  number  as  a  function  of  distance 
from  the  wall  located  at  40  cm  for  Mach  15.9  case. 


A.  Fluctuation  Patterns 

The  CR  model  depends  on  the  values  of  the  cross- 
sections  of  elementary  processes.  Some  of  them  are  un¬ 
certain,  in  particular  atom-atom  inelastic  collisions;  thus, 
in  article  1  the  CR  model  was  calibrated  against  experi¬ 
mental  data  for  the  steady-state  case.  In  this  section  we 
examine  the  effect  of  unsteady  coupling  between  the  ki¬ 
netics  and  transport,  using  the  unsteady  solver  described 
above.  Initial  work  on  such  unsteady  effects  verified  the 
exsistence  of  oscillations  in  the  form  of  longitudinal  pul¬ 
sations  in  the  direction  of  the  flow5.  Building  upon  this 
work,  similar  numerical  experiments  are  carried  out  here, 
albeit  with  a  more  detailed  CR  model  and  higher-order 
accuracy. 

The  conditions  of  the  shock  tube  experiments  were 
simulated  numerically  by  impinging  argon  gas  upon  an 
ideally  reflecting  wall,  thereby  initiating  a  shock  that 
propagates  against  the  incoming  flow  as  illustrated  in 
Figure  1.  Note  that  the  free-stream  conditions  listed 
therein  correspond  to  a  shock  with  an  initial  strength  of 
Mach  19.  As  the  shock  reflects  from  the  wall,  its  Mach 
number  decreases  as  the  thermal  energy  is  quickly  con¬ 
verted  to  electronic  energy  in  the  form  of  bound  excited 
states  as  well  as  free  electrons,  until  a  quasi-steady  limit 
is  reached.  Thus,  the  freestream  conditions  provided  in 
the  figure  are  those  used  to  reproduce  the  Mach  15.9 
shock  in  the  laboratory  frame,  corresponding  to  the  first 
case  of  Table  I  of  Part  1.  The  computational  domain  be¬ 
ing  fixed,  a  higher  resolution  can  be  obtained  for  choos¬ 
ing  this  reference  frame  with  a  lower  propagation  speed 


of  the  shock  front,  while  avoiding  the  potential  effects  of 
a  subsonic  boundary  condition.  During  the  calculations, 
we  can  easily  monitor  local  and  global  variables,  such  as 
the  instantaneous  shock  Mach  number  in  the  laboratory 
frame  and  the  ”  induction”  length,  i.e.  the  separation  be¬ 
tween  the  shock  front  and  the  avalanche  region;  the  latter 
is  identified  by  the  location  where  the  electron  density  has 
reached  peak  value. 

The  results  in  Figures  2  and  3  clearly  show  undamped, 
periodic  fluctuations  in  both  the  shock  Mach  number  as 
well  as  the  induction  length  (in  blue)  compared  to  the 
experimental  (averaged)  value  (dashed  red  line).  These 
oscillations  have  a  periodicity  of  approximately  32.5  //sec 
and  are  characteristic  of  the  experimental  observations. 
As  evident  in  the  discontinuous  shifts  of  the  induction 
length,  the  oscillations  do  not  vary  smoothly,  signal¬ 
ing  a  strong  non-linear  behavior.  In  fact,  this  behavior 


FIG.  3.  Instantaneous  induction  length  as  a  function  of  time 
for  Mach  15.9  case.  Dashed  red  line  indicates  the  experimen¬ 
tally  observed  value  of  2  cm.  The  fluctuation  periodicity  is 
32.5  //sec,  with  an  induction  length  of  2.22  cm  averaged  over 
one  period. 
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FIG.  4.  Dynamic  evolution  of  ID  ionizing  shock  structure:  a) 
quasi-equilibrium  state  b)  pressure  wave  initiated  at  electron 
avalanche  travels  towards  compression  shock  c)  after  overtak¬ 
ing  shock,  pressure  wave  is  reflected  as  an  entropy  wave  due  to 
strengthening  of  shock  d)  sensitivity  of  excitation  and  ioniza¬ 
tion  rates  to  temperature  jump  across  entropy  wave  results  in 
earlier  onset  of  electron  avalanche.  The  process  then  repeats. 

could  be  well  approximated  by  the  square-wave  model 
for  detonations10  in  which  the  kinetics  are  assumed  to 
be  frozen  below  a  certain  threshold,  but  infinitely  fast 
above  it,  reaching  equilibrium  instantaneously.  This  ne¬ 
cessitates  the  introduction  of  an  average  induction  length 
over  a  period, 


with  being  the  instantaneous  induction  length.  The 
computed  value  is  shown  in  3  and  compares  well  with 
the  experimental  value  (itself  an  approximation). 

In  regards  to  the  discontinuous  change  in  the  induction 
length  exhibited  in  this  particular  case,  there  is  evidence 
that  these  anomalies  have  a  physical  origin.  A  sequence 
of  characteristic  profiles  is  shown  in  Figure  4.  From  top- 
down  and  left-right,  the  computed  profiles  of  the  mass 
density  (red  line)  and  electron  density  (blue  line)  are 
shown,  compared  to  a  snapshot  of  the  experimental  val¬ 
ues  (symbols).  Disturbances  of  the  total  and  electron 
densities  are  seen  to  propagate  from  the  avalanche  region 
towards  the  shock  and  reflected  back.  The  perturbation 
in  temperature  is  such  that,  given  the  highly  non-linear 
kinetics  within  the  induction  zone  (see  article  1),  the  lo¬ 
cation  of  the  avalanche  can  be  significantly  altered.  As 
seen  in  Figure  4-d,  a  new  ’’spot”  where  the  ionization 


avalanche  occurs  can  be  clearly  identified.  This  will  cre¬ 
ate  in  turn  a  new  location  for  the  peak  ionization  and 
a  corresponding  pressure  wave  (electron  pressure)  which 
couples  back  to  the  shock  front,  and  the  process  is  re¬ 
peated. 

The  nature  of  this  process  can  be  better  demonstrated 
with  the  space-time  plots  given  in  Figure  5.  The  plots 
clearly  show  pressure  waves  emerging  from  the  electron 
avalanche,  that  travel  in  both  directions.  This  can  be 
considered  the  initiation  of  the  fluctuation  cycle.  The 
forward-moving  wave  (left  in  the  x-t  diagram)  overtakes 
the  shock,  thereby  increasing  its  temperature  and  speed; 
the  interaction  also  produces  a  rarefaction  wave  due  to 
the  acceleration  of  the  shock  which  is  too  weak  to  notice 
on  that  figure.  The  resulting  high  temperature  region 
generated  behind  the  shock  is  separated  from  the  previ¬ 
ous  post-shock  region  by  a  contact  wave  that  can  be  seen 
as  traveling  away  from  the  shock  as  it  accelerates.  Identi¬ 
fication  of  the  contact  wave  is  trivially  done  by  compar¬ 
ing  the  x-t  diagrams  for  pressure  and  density.  A  more 
detailed  reproduction  of  the  (computed)  x-t  diagram  is 
shown  in  Figure  6. 

The  temperature  jump  across  the  contact  wave  is  cru¬ 
cial.  Due  to  the  exponential  dependence  on  temperature, 
a  relatively  small  jump  can  accelerate  the  kinetics  and 
cause  the  electron  avalanche  to  form  much  sooner  than 
its  quasi-equilibrium  position.  If  the  non-linear  process 
leading  to  the  avalanche  is  sufficiently  fast  and  sensitive, 
this  shift  in  the  avalanche  region  can  be  discontinuous, 
as  seen  in  Figure  4,  i.e.  the  electron  density  profile  in  the 
induction  region  is  non-monotonous — alternatively,  if  we 
consider  the  location  of  peak  ionization,  we  see  a  bifurca¬ 
tion.  This  process  is  entirely  analogous  to  the  mechanism 
for  detonation  wave  fluctuations  proposed  by  Alpert  and 
Toong11,  and  was  reproduced  in  numerical  simulations  of 
gaseous  detonations  (e.g.  see  12). 

This  explanation  of  the  fluctuation  mechanism  can  be 
confirmed  by  a  simple  estimation  of  the  periodicity  of 
the  oscillations  from  basic  wave  theory.  In5  it  was  shown 
that  the  period  can  be  approximated  by 

r  =  l(—^—  +  —),  (16) 

\a2-u2  u2J 

assuming  that  the  pressure  wave  travels  towards  the 
shock  with  the  nonlinear  wave  speed  a2  —  u2,  while  the 
entropy  wave  reflected  from  the  shock  travels  with  ve¬ 
locity  u2 ,  where  the  subscript  2  denotes  the  post-shock 
state.  With  u2  =  1284m/s  and  a2  =  2852m/s  taken  as 
the  post-shock  fluid  and  sound  speeds  respectively,  the 
estimated  periodicity  is  31.4//sec  which  agrees  well  with 
the  observed  value  of  32.5/isec  (cf.  Figure  3). 


B.  Sensitivity  to  cross  sections 

We  examined  the  calibration  of  some  key  cross-sections 
in  article  1  using  the  induction  length  in  the  steady-state 
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FIG.  5.  x-t  diagrams  for  ID  Ma  14.7  case.  From  top  row: 
total  density  (p)  on  the  left,  total  pressure  (p)  on  the  right; 
heavy-particle  temperature  (T^,  left),  electron  temperature 
(Te,  right);  on  the  last  row,  velocity  (14,  left)  and  ionization 
fraction  (a,  right).  Gradients  are  enhanced  through  shading 
effects.  Colormap:  min  H  I  max 


FIG.  6.  x-t  diagrams  for  ID  Ma  14.7  case  -  total  density 
(detail) . 
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FIG.  7.  Sensitivity  of  induction  length  and  shock  Mach 
number  to  cross  sections  for  ID  unsteady  calculations;  using 
Drawin’s  excitation  cross  sections. 


FIG.  8.  Sensitivity  of  induction  length  and  shock  Mach  num¬ 
ber  to  cross  sections  for  ID  unsteady  calculations;  reducing 
the  atom  impact  excitation  cross-sections  by  factor  of  10 


case,  but  it  may  also  be  possible  to  use  the  unsteady  fea¬ 
tures  to  place  bounds  on  these  cross-sections.  For  that, 
knowing  the  effect  of  the  cross  sections  on  the  dynamics  of 
the  oscillations  is  crucial.  While  an  analytical  approach 
to  stability  analysis  is  not  pursued  here,  semi-empirical 
considerations  can  serve  as  a  valuable  precursor  to  such 
studies. 

As  a  first  test,  the  Mach  14.7  case  was  run  using  the 
atom-atom  impact  excitation  cross  sections  as  computed 
from  Drawin’s  formula,  Eq.  (17),  without  modification. 


<  (£)  =  4™o  (  — 


mH 


<2f ; 


2  me 


e/Ii  -  1 


rriAr  +  me 


1  + 


2m  e 


mAr+rnt 


ie/Ii-  1))' 


(17) 


where  Ih  is  the  ionization  potential  of  the  hydrogen  atom 
and  £  =  6  is  the  number  of  optical  electrons  of  argon.  As 
a  reminder,  the  resulting  rates  are  more  than  a  factor  of 
10  greater  than  those  used  to  obtain  the  solutions  with 
the  experimentally-observed  induction  length.  Compar¬ 
ing  the  results  shown  in  Figures  7  and  8  indicate  that 
the  large  atom-impact  excitation  cross  sections  obtained 
from  Drawin’s  formula  essentially  halves  the  induction 
length  compared  to  the  observed  value  (8).  Furthermore, 
the  influence  on  the  oscillations  is  quite  pronounced,  re¬ 
sulting  in  significant  damping  to  the  point  that  the  oscil¬ 
lations  can  no  longer  be  sustained.  This  can  be  explained 
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as  follows:  as  the  pressure  wave  emanating  from  the 
electron  avalanche  overtakes  the  compression  shock,  the 
heavy  particle  temperature  increases.  This  in  turn  in¬ 
creases  the  atom-impact  excitation,  the  rate  having  a  sen¬ 
sitive  exponential  dependence  on  T^.  The  atom-impact 
inelastic  processes  acts  as  an  energy  sink,  effectively  ab¬ 
sorbing  energy  from  the  pressure  wave  and  dampening 
the  oscillations. 

In  a  separate  test,  the  electron- atom  impact  excitation 
cross  sections  were  reduced  by  a  factor  of  10.  In  this 
case  the  atom  impact  excitation  cross  sections  were  also 
modified  so  that  the  induction  length  was  close  to  4.4  cm 
as  seen  in  Figure  8.  While  the  oscillations  still  persisted 
in  this  case,  the  discontinuous  change  in  the  induction 
length  (cf.  Figure  2)  had  disappeared.  Instead,  the  oscil¬ 
lations  vary  smoothly  over  the  entire  period  of  the  fluc¬ 
tuation.  This  key  observation  can  be  explained  from  the 
association  of  the  discontinuous  shift  with  a  highly  sen¬ 
sitive  non-linear  coupling  with  exponential  temperature 
dependence  of  the  rate  equations;  this  dependence  can 
be  well- approximated  by  the  square- wave  model  for  det¬ 
onation  phenomena  as  previously  noted.  As  the  electron- 
impact  cross  sections  are  decreased,  the  strength  of  the 
non-linear  coupling  diminishes  and  other  processes  such 
as  thermalization  become  relatively  stronger.  This  slows 
the  rate  of  onset  of  the  electron  avalanche  and  gives  the 
fluid  more  time  to  “react”  to  the  kinetics. 


IV.  2D  RESULTS 
A.  Test  Setup 

With  a  general  understanding  of  the  oscillation  mech¬ 
anism  and  the  resulting  wave  pattern  obtained  from  the 
unsteady  ID  simulations,  we  can  now  extend  the  simula¬ 
tions  and  theory  to  a  two-dimensional  problem.  With  an¬ 
other  dimension,  the  observed  oscillations  are  no  longer 
confined  to  the  longitudinal  direction,  leading  to  the  pos¬ 
sibility  of  transverse  waves.  While  this  opens  the  door 
to  more  complex  wave  interactions  and  structures,  a  sys¬ 
tematic  approach  that  extends  the  ideas  of  the  previous 
section  in  a  straight-forward  manner  is  taken  here. 

The  2D  simulations  were  conducted  on  the  domain 
0  <  x  <  54  cm,  0  <  y  <  18  cm,  with  a  depth  of 
10  cm13.  A  general  schematic  is  provided  in  Figure  9 
and  the  simulations  were  initiated  in  the  same  manner 
as  the  ID  simulations,  with  a  high  speed  flow  imping¬ 
ing  upon  a  wall,  resulting  in  a  shock  propagating  in  the 
opposite  direction.  The  difference  being,  however,  that 
the  2D  simulations  required  an  initial  perturbation  to 
initiate  disturbances  in  the  transverse  direction.  While 
such  disturbances  can  have  any  number  of  origins  in  an 
actual  shock  tube14,  they  are  explicitly  included  in  the 
numerical  experiment  as  small  perturbations  in  order  to 
accelerate  the  growth  of  instabilities  and  make  the  nu¬ 
merical  simulations  tractable.  Such  perturbations  must 
be  sufficiently  random  and  weak  as  to  not  overshadow 


54  cm  10  cm 

V 


FIG.  9.  Numerical  setup  for  2D  unsteady  shock  tube  simula¬ 
tions. 

the  dynamics  inherent  in  the  model,  yet  strong  enough 
to  initiate  pressure  disturbances  that  can  essentially  be 
amplified  by  the  true  dynamics. 

Simulations  have  been  performed  for  two  of  the  cases, 
Mach  16.5  and  14.7  .  All  simulations  shown  were  run 
on  a  Cartesian  mesh  with  Ax  =  Ay  =  0.33  mm  and 
have  taken  into  account  excited  levels  only  from  the  4 s 
manifold  of  neutral  argon.  Although  it  was  shown  in 
the  ID  simulations  that  higher  levels  are  required  to  ob¬ 
tain  good  agreement  in  the  radiative  cooling  region,  the 
induction  zone  was  for  the  most  part  unaffected.  There¬ 
fore,  neglecting  the  levels  beyond  those  contained  in  the 
4 s  manifold  should  have  negligible  effect  on  the  overall 
dynamics  of  the  instability  in  the  induction  zone,  which  is 
the  primary  interest  here,  while  significantly  accelerating 
the  computations. 

B.  Mach  16.5  case 

Figure  10  shows  the  evolution  of  the  Mach  16.5  ioniz¬ 
ing  shock  from  a  nearly  planar  wave  to  a  self-sustaining 
oscillating  pattern.  Contours  of  the  computed  refractive 
index  are  shown,  expecially  useful  to  visualize  the  density 
gradients  with  a  high  level  of  contrast.  The  initial  pertur¬ 
bation,  visible  as  a  half-sine  wave  bulge  in  the  shock  front, 
generates  transverse  waves  that  are  apparently  random 
in  nature.  After  propagating  a  few  centimeters,  however, 
the  waves  begin  to  oscillate  with  a  resonant  frequency, 
creating  the  strong  incident  and  reflected  shock  pattern 
that  is  clearly  visible. 

Figure  11  provides  snapshots  of  the  same  simulation 
at  a  later  time  as  generated  by  two  different  visualiza¬ 
tion  techniques. The  top  figure  shows  a  finite-fringe  inter- 
ferogram  generated  from  the  numerical  simulations;  this 
clearly  shows  the  corrugations  in  the  shock  front  and  an 
overall  pattern  similar  to  the  UTIAS  experiments.  The 
refractive  index  plot  below,  taken  at  a  later  time,  fur¬ 
ther  reveals  a  well-defined  structure  behind  the  shock, 
having  a  periodic  nature  in  the  transverse  direction,  and 
strong  vorticity  generated  within  the  shock  region.  From 
an  investigation  of  the  wave  structure  it  was  determined 
that  the  shock  front  actually  consists  of  incident  and  re¬ 
flected  shocks  as  well  as  Mach  stems,  the  intersection  of 
which  forms  a  series  of  triple  points.  These  nodes  appear 
in  a  quite  regular  pattern  in  both  the  longitudinal  and 
transverse  directions,  indicating  a  resonant  phenomenon. 
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FIG.  10.  Evolution  of  shock  structure  for  Ma  16.5.  Contours 
of  refractive  index. 


FIG.  11.  Simulated  interferogram  and  refractive  index  for  2D 
Mach  16.5  shock.  The  shock  front  consists  of  incident  and 
reflected  shock  waves.  Color  map:  min  H  M  max 


FIG.  12.  Trace  of  triple  points  for  Mach  16.5  case  expos¬ 
ing  ionizing  cell  structure.  Brighter  regions  correspond  to 
stronger  flow  vorticity. 


FIG.  13.  Simulated  interferogram  and  refractive  index  for  2D 
Mach  14.7  shock.  Colormap:  min  H  1  max 

These  triple  points  are  the  source  of  the  vorticity  ob¬ 
served  as  far  downstream  as  10-15  times  the  induction 
length.  The  trace  of  these  triple  points,  given  in  Figure 
12,  highlights  a  cellular  structure  remarkably  similar  to 
detonation  cells.  These  analogous  ionization  cells  initiate 
at  the  intersection  of  two  triple  points,  which  are  not  only 
the  source  of  maximum  vorticity,  but  also  of  maximum 
temperature  and  ionization;  they  can  be  thought  as  the 
two-dimensional  analog  to  the  one-dimensional  electron 
avalanche.  The  resulting  pressure  wave  generated  by  the 
increase  in  electron  number  density  at  this  point  trav¬ 
els  radially  outward,  as  indicated  by  the  trace  of  triple 
points.  Figure  12  can  be  considered  the  equivalent  of 
a  soot  trace  commonly  used  in  experimental  studies  of 
gaseous  detonations. 


C.  Mach  14.7  case 

The  case  of  a  Mach  14.7  shock  was  also  simulated,  and 
from  the  snapshots  in  Figure  13,  as  well  as  the  simulated 
soot  trace  in  Figure  14,  it  is  clear  that  decreasing  the 
Mach  number  has  increased  the  cell  size  in  both  the  lon¬ 
gitudinal  and  transverse  directions.  A  direct  comparison 
of  the  experimental  and  simulated  results  for  the  Mach 
14.7  case  is  presented  in  Figure  15.  As  it  can  be  seen, 
the  shock  structure  is  well-predicted  by  the  numerical 


FIG.  14.  Numerical  soot  trace  for  Mach  14.7  case,  exposing 
ionizing  cell  structure.  Image  was  made  by  super-imposing 
instantaneous  snapshots  of  the  shock  front. 


FIG.  15.  Direct  comparison  of  experimental  and  numerical 
results  for  Mach  14.7  shock.  Clearly  evident  are  the  corru¬ 
gations  in  the  compression  shock  as  well  as  unique  formation 
patterns  in  the  electron  avalanche  which  are  accurately  repro¬ 
duced  by  the  solver. 


simulations,  with  similar  relaxation  lengths  and  discon¬ 
tinuous  gradients  in  the  electron  avalanche  region  along 
the  transverse  direction.  The  agreement  is  sufficiently 
close  that  a  similar  shock  structure  can  be  overlaid  upon 
on  the  experimental  results  showing  excellent  agreement 
and  compatibility  with  the  visible  structures.  In  partic¬ 
ular,  the  pattern  of  incident  and  reflected  shocks  is  clear 
in  the  fringe  shifts  of  the  interferogram. 


D.  Cellular  structure 

Estimation  of  the  2D  cell  size  was  provided  in5,  us¬ 
ing  the  ID  periodicity  calculation  and  assuming  that  the 
longitudinal  periodicity  must  equal  the  transverse  wave 
propagation  time  across  the  height  (y- dir)  of  one  com¬ 
plete  cell.  The  cell  width  could  therefore  be  estimated 
by 


5  ~  (a2  +  p-fjr,  (18) 

where  the  bulk  velocity  in  the  transverse  direction  behind 
the  shock,  /;2,  is  assumed  to  be  zero,  and  S  corresponds  to 
a  half-wavelength  of  the  shock  curvature  observed  in  the 
UTIAS  experiments.  Using  the  periodicity  as  determined 


FIG.  16.  Details  of  the  2D  ionizing  shock  structure  for  the 
Mach  14.7  case. 


in  the  ID  Mach  15.9  simulations  as  a  reference  case,  the 
periodicity  of  the  Mach  14.7  and  16.5  cases  can  roughly 
be  approximated  by 

,  Ma  t' 

T  ~ 

From  this  it  is  found  r'Mal47  ~  77  //sec  and  r'Mal6  5  ~ 
28  // sec.  Inserting  these  values  into  Eq.  18  gives 
^Mai4.7  —  (2658  m/s) (77  //sec)  ~  20  cm  and  5mcl16.5  — 
(2973  m/s)  (28  //sec)  ~  8.3  cm.  The  predicted  cell 
height  for  the  Mach  16.5  case  agrees  favorably  with  the 
numerically-observed  value  of  7.2  cm,  while  the  Mach 
14.7  case  is  somewhat  over-predicted,  with  a  maximum 
observed  value  of  just  under  14  cm.  However,  this  case  is 
the  most  sensitive  to  finite-size  effects  of  the  shock  tube; 
as  we  expect  boundary-layer  effects  to  shorten  the  effec¬ 
tive  induction  length,  the  transverse  size  of  the  ioniza¬ 
tion  cell  would  be  also  affected.  Notwithstanding  these 
boundary-layer  effects,  we  can  claim  a  good  agreement 
over  the  range  of  Mach  number  studied. 

The  experimental  shock  structure  can  be  inferred  from 
the  simulated  interferogram.  In  IV  D,  we  trace  the  inci¬ 
dent  and  reflected  shocks  matching  the  shift  in  fringes  in 
the  induction  zone. 

In  a  similar  fashion,  Buchanan15  drew  schematics  of 
the  Schlieren  images,  reproduced  here  in  Fig.  17.  These 
also  indicate  the  ionization  cell  structure  and  its  size.  It 
is  of  course  more  difficult  to  manually  extract  this  in¬ 
formation  from  Schlieren  pictures  than  from  numerical 
simulations.  Nevertheless,  Buchanan  was  able  to  deduce 
a  trend  for  the  correlation  between  the  number  of  ob¬ 
served  half- wavelengths  (cell  size),  the  equilibrium  val¬ 
ues  of  the  electron  number  density  (n*)  and  the  shock 
Mach  number  ( Ms ).  Buchanan  postulated  that  the  en¬ 
ergy  transfer  rate  between  electrons  and  heavy  particles 
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(a)Ma  =  12.71, 
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(c)Ma  =  14.38, 
pi  =  4.11,  T\  =  297.0 


(d)Ma  =  16.10,  pi  =  5.79, 
T\  =  299.4 


FIG.  17.  Schematics  of  experimental  Schlieren  images  ob¬ 
tained  by  Buchanan15.  The  shock  structure  to  the  right  of 
the  translational  shock  front  (thick  dotted  line)  clearly  show 
the  structure  of  the  ionization  cells.  The  number  of  half-cells 
per  shocktube  channel  width  is  2,  5,  3,  and  6  for  cases  (a)-(d) 
respectively. 


played  a  role  in  the  oscillations,  citing  the  work  of  Morse 
and  Ingard16,  and  used  —  Te,  ne,  n h  as  indicators.  As 
was  shown  in  Part  1  with  steady-state  calculations,  this 
role  is  played  by  the  ladder-climbing  mechanism,  and  Te 
is  only  an  approximate  indication  of  the  status  of  the 
Atomic  State  Distribution  Function  (ASDF),  which  is 
characterized  by  two  sets  of  temperatures,  describing  the 
deviations  from  Boltzmann  and  Saha  equilibrium.  Thus, 
the  picture  is  somewhat  more  complex  than  Buchanan 
envisioned.  Nevertheless,  it  would  be  worthwhile  at  some 
point  in  the  future  to  pursue  a  systematic  study  of  all 
cases  considered  by  Buchanan  and  compare  the  results. 


V.  CONCLUSIONS 


Following  Part  1  of  this  article  where  a  detailed 
collisional-radiative  model  of  argon  was  described  and 
tested  against  structures  of  ionizing  shocks  in  argon, 
we  extended  the  computations  to  account  for  unsteady 
and  multi-dimensional  effects.  We  found  that  the  struc¬ 
ture  of  density  profiles  within  the  induction  zone  can 
be  very  easily  explained  by  a  periodic  wave  propaga¬ 
tion  between  the  avalanche  region  and  the  shock  front. 
Fluctuations  in  electron  pressure  (and  therefore  total 
pressure)  emanate  from  the  avalanche  region,  subject 
to  very  rapid  changes  in  electron  kinetics  and  propa¬ 
gate  towards  the  shock,  where  they  are  reflected.  The 
shock  is  accelerated  and  the  reflected  contact  wave  sepa¬ 
rates  a  region  of  higher  temperature  and  faster  kinetics, 
leading  to  a  new  fluctuation  in  the  location  of  the  ion¬ 
ization  avalanche.  When  mutli-dimensional  effects  are 
considered,  this  process  invites  the  formation  of  trans¬ 
verse  waves,  and  a  cellular  structure  is  clearly  observed. 
This  pattern  matches  experimental  observations  fairly 
well,  accounting  for  limitations  in  both  the  numerical 
model  (no  boundary  layer,  2D  only)  and  the  experi¬ 
ments.  This  formerly  conjectured5  dynamical  effect  has 
now  been  convincingly  demonstrated  here,  as  well  as  the 
analogy  with  cellular  structures  in  gaseous  detonations, 
as  evidenced  for  example  by  the  numerical  soot  trace 
patterns.  Thus,  the  observation  of  such  structures  in 
non-equilibrium  shock  dynamics  does  not  have  to  rely  on 
the  exo-thermicity  of  chemical  reactions,  but  on  the  pres¬ 
ence  of  a  significant  induction  zone  followed  by  a  narrow 
avalanche  region  characterized  by  very  rapid,  non-linear 
chemical  kinetics,  and  thus  extremely  sensitive  to  fluctu¬ 
ations  in  the  conditions  within  the  induction  zone.  The 
ladder-climbing  process  of  a  noble  gas  ASDF  provides 
such  conditions. 

There  is  clearly  a  lot  that  can  be  learned  from  shock- 
tube  studies,  and  we  have  shown  that  for  argon,  the  un¬ 
steady  phenomena,  if  properly  visualized,  can  provide 
a  wealth  of  information  to  better  validate  the  chemical 
(collisional-radiative)  models.  This  also  implies  that  the 
corresponding  numerical  models  must  be  able  to  accu¬ 
rately  reproduce  these  unsteady  and  multi-dimensional 
effects.  As  shown  here,  this  capability  is  within  reach. 
Although  we  have  limited  the  level  of  detail  of  the  CR 
model  for  the  multi-dimensional  studies,  we  have  a  strong 
confidence  that  the  detailed  model  can  be  applied  in  2D 
and  even  3D  simulations  on  reasonably-sized  computing 
platforms.  We  are  currently  working  towards  improv¬ 
ing  the  efficiency  of  the  computational  approach,  and 
since  most  of  the  computational  cost  lies  in  the  CR  ki¬ 
netics,  higher  dimensionality  or  addition  of  the  Navier- 
Stokes  terms  does  not  present  undue  challenges.  Thus, 
future  work  will  aim  at  more  complete  and  more  detailed 
capabilities  in  those  directions;  combined  with  state-of- 
the-art  visualization  and  diagnostic  techniques  with  high 
time  and  spatial  resolution,  such  numerical  capabilities 
will  prove  valuable  for  the  study  of  non-equilibrium  high- 


10 


temperature  gas  dynamics  in  modern  shock-tube  exper¬ 
iments. 
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where  s  is  the  species  index,  ms  its  molecular  mass  and 
cv,s  the  constant- volume  specific  heat  (per  unit  mass)  of 
that  species.  The  total  pressure  derivatives,  which  enter 
in  the  definition  of  the  flux  jacobian,  are: 


pe=(^)^=7',_i  <a2) 
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